{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">Registration: Memory-Time Trade-off</h1>\n",
    "\n",
    "When developing a registration algorithm or when selecting parameter value settings for an existing algorithm our choices are dictated by two, often opposing, constraints:\n",
    "<ul>\n",
    "<li>Required accuracy.</li>\n",
    "<li>Alloted time.</li>\n",
    "</ul>\n",
    "\n",
    "As the goal of registration is to align multiple data elements into the same coordinate system, it is only natural that the primary focus is on accuracy. In most cases the reported accuracy is obtained without constraining the algorithm's execution time. Don't forget to provide the running times even if they are not critical for your particular application as they may be critical for others. \n",
    "\n",
    "With regard to the emphasis on execution time, on one end of the spectrum we have longitudinal studies where time constraints are relatively loose. In this setting a registration taking an hour may be perfectly acceptable. At the other end of the spectrum we have intra-operative registration. In this setting, registration is expected to complete within seconds or minutes. The  underlying reasons for the tight timing constraints in this setting have to do with the detrimental effects of prolonged anesthesia and with the increased costs of operating room time. While short execution times are important, simply completing the registration on time without sufficient accuracy is also unacceptable.  \n",
    "\n",
    "This notebook illustrates a straightforward approach for reducing the computational complexity of registration for intra-operative use via preprocessing and increased memory usage, a case of the [memory-time trade-off](https://en.wikipedia.org/wiki/Space%E2%80%93time_tradeoff). \n",
    "\n",
    "The computational cost of registration is primarily associated with interpolation, required for evaluating the similarity metric. Ideally we would like to use the fastest possible interpolation method, nearest neighbor. Unfortunately, nearest neighbor interpolation most often yields sub-optimal results. A straightforward solution is to pre-operatively create a super-sampled version of the moving-image using higher order interpolation*. We then perform registration using the super-sampled image, with nearest neighbor interpolation.\n",
    "\n",
    "Tallying up time and memory usage we see that:\n",
    "\n",
    "<table>\n",
    "  <tr><td></td> <td><b>time</b></td><td><b>memory</b></td></tr>\n",
    "  <tr><td><b>pre-operative</b></td> <td>increase</td><td>increase</td></tr>\n",
    "  <tr><td><b>intra-operative</b></td> <td>decrease</td><td>increase</td></tr>\n",
    "</table><br><br>  \n",
    "\n",
    "\n",
    "<font size=\"-1\">*A better approach is to use single image super resolution techniques such as the one described in A. Rueda, N. Malpica, E. Romero,\"Single-image super-resolution of brain MR images using overcomplete dictionaries\", <i>Med Image Anal.</i>, 17(1):113-132, 2013.</font> \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "import SimpleITK as sitk\n",
    "\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#utility method that either downloads data from the Girder repository or\n",
    "#if already downloaded returns the file name for reading from disk (cached data)\n",
    "%run update_path_to_download_script\n",
    "from downloaddata import fetch_data as fdata\n",
    "\n",
    "import registration_utilities as ru\n",
    "\n",
    "from ipywidgets import interact, fixed\n",
    "\n",
    "\n",
    "def register_images(fixed_image, moving_image, initial_transform, interpolator):\n",
    "\n",
    "    registration_method = sitk.ImageRegistrationMethod()    \n",
    "    registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)\n",
    "    registration_method.SetMetricSamplingStrategy(registration_method.REGULAR)\n",
    "    registration_method.SetMetricSamplingPercentage(0.01)\n",
    "    registration_method.SetInterpolator(interpolator)     \n",
    "    registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=1000) \n",
    "    registration_method.SetOptimizerScalesFromPhysicalShift() \n",
    "    registration_method.SetInitialTransform(initial_transform, inPlace=False)\n",
    "    \n",
    "    final_transform = registration_method.Execute(fixed_image, moving_image)\n",
    "    \n",
    "    return( final_transform, registration_method.GetOptimizerStopConditionDescription())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data\n",
    "\n",
    "We use the the training data from the Retrospective Image Registration Evaluation (<a href=\"http://www.insight-journal.org/rire/\">RIRE</a>) project.\n",
    "\n",
    "The RIRE reference, ground truth, data consists of a set of corresponding points in the fixed and moving coordinate systems. These points were obtained from fiducials embedded in the patient's skull and are thus sparse (eight points). We use these to compute the rigid transformation between the two coordinate systems, and then generate a dense reference. This generated reference data is more similar to the data you would use for registration evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fixed_image =  sitk.ReadImage(fdata(\"training_001_ct.mha\"), sitk.sitkFloat32)\n",
    "moving_image = sitk.ReadImage(fdata(\"training_001_mr_T1.mha\"), sitk.sitkFloat32) \n",
    "fixed_fiducial_points, moving_fiducial_points = ru.load_RIRE_ground_truth(fdata(\"ct_T1.standard\"))\n",
    "\n",
    "R, t = ru.absolute_orientation_m(fixed_fiducial_points, moving_fiducial_points)\n",
    "reference_transform = sitk.Euler3DTransform()\n",
    "reference_transform.SetMatrix(R.flatten())\n",
    "reference_transform.SetTranslation(t)\n",
    "\n",
    "# Generate a reference dataset from the reference transformation (corresponding points in the fixed and moving images).\n",
    "fixed_points = ru.generate_random_pointset(image=fixed_image, num_points=1000)\n",
    "moving_points = [reference_transform.TransformPoint(p) for p in fixed_points]    \n",
    "\n",
    "interact(lambda image1_z, image2_z, image1, image2:ru.display_scalar_images(image1_z, image2_z, image1, image2), \n",
    "         image1_z=(0,fixed_image.GetSize()[2]-1), \n",
    "         image2_z=(0,moving_image.GetSize()[2]-1), \n",
    "         image1 = fixed(fixed_image), \n",
    "         image2=fixed(moving_image));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Invest time and memory in exchange for future time savings\n",
    "\n",
    "We now resample our moving image to a finer spatial resolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Isotropic voxels with 1mm spacing.\n",
    "new_spacing = [1.0]*moving_image.GetDimension()\n",
    "\n",
    "# Create resampled image using new spacing and size.\n",
    "original_size = moving_image.GetSize()\n",
    "original_spacing = moving_image.GetSpacing()\n",
    "resampled_image_size = [int(spacing/new_s*size) \n",
    "                        for spacing, size, new_s in zip(original_spacing, original_size, new_spacing)]  \n",
    "resampled_moving_image = sitk.Image(resampled_image_size, moving_image.GetPixelID())\n",
    "resampled_moving_image.SetSpacing(new_spacing)\n",
    "resampled_moving_image.SetOrigin(moving_image.GetOrigin())\n",
    "resampled_moving_image.SetDirection(moving_image.GetDirection())\n",
    "\n",
    "# Resample original image using identity transform and the BSpline interpolator.\n",
    "resample = sitk.ResampleImageFilter()\n",
    "resample.SetReferenceImage(resampled_moving_image)                      \n",
    "resample.SetInterpolator(sitk.sitkBSpline)  \n",
    "resample.SetTransform(sitk.Transform())\n",
    "resampled_moving_image = resample.Execute(moving_image)\n",
    "\n",
    "print('Original image size and spacing: {0} {1}'.format(original_size, original_spacing)) \n",
    "print('Resampled image size and spacing: {0} {1}'.format(resampled_moving_image.GetSize(), \n",
    "                                                         resampled_moving_image.GetSpacing()))\n",
    "print('Memory ratio: 1 : {0}'.format((np.array(resampled_image_size)/np.array(original_size).astype(float)).prod())) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another option for resampling an image, without any transformation, is to use the ExpandImageFilter or \n",
    "in its functional form SimpleITK::Expand. This filter accepts the interpolation method and an integral expansion factor. This is less flexible than the resample filter as we have less control over the resulting image's spacing. \n",
    "On the other hand this requires less effort from the developer, a single line of code as compared to the cell above:\n",
    "\n",
    "resampled_moving_image = sitk.Expand(moving_image, \n",
    "                                     [int(original_s/new_s + 0.5) for original_s, new_s in zip(original_spacing, new_spacing)], \n",
    "                                     sitk.sitkBSpline)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Registration\n",
    "\n",
    "### Initial Alignment\n",
    "\n",
    "We will use the same initial alignment for both registrations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "initial_transform = sitk.CenteredTransformInitializer(fixed_image, \n",
    "                                                      moving_image, \n",
    "                                                      sitk.Euler3DTransform(), \n",
    "                                                      sitk.CenteredTransformInitializerFilter.GEOMETRY)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Original Resolution\n",
    "\n",
    "For this registration we use the original resolution and linear interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%timeit -r1 -n1\n",
    "# The arguments to the timeit magic specify that this cell should only be run once. \n",
    "\n",
    "# We define this variable as global so that it is accessible outside of the cell (timeit wraps the code in the cell\n",
    "# making all variables local, unless explicitly declared global).\n",
    "global original_resolution_errors\n",
    "\n",
    "final_transform, optimizer_termination = register_images(fixed_image, moving_image, initial_transform, sitk.sitkLinear)\n",
    "final_errors_mean, final_errors_std, _, final_errors_max, original_resolution_errors = ru.registration_errors(final_transform, fixed_points, moving_points)\n",
    "\n",
    "print(optimizer_termination)\n",
    "print('After registration, errors in millimeters, mean(std): {:.2f}({:.2f}), max: {:.2f}'.format(final_errors_mean, final_errors_std, final_errors_max))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Higher Resolution\n",
    "\n",
    "For this registration we use the higher resolution image and nearest neighbor interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%timeit -r1 -n1\n",
    "# The arguments to the timeit magic specify that this cell should only be run once. \n",
    "\n",
    "# We define this variable as global so that it is accessible outside of the cell (timeit wraps the code in the cell\n",
    "# making all variables local, unless explicitly declared global).\n",
    "global resampled_resolution_errors \n",
    "\n",
    "final_transform, optimizer_termination = register_images(fixed_image, resampled_moving_image, initial_transform, sitk.sitkNearestNeighbor)\n",
    "final_errors_mean, final_errors_std, _, final_errors_max, resampled_resolution_errors = ru.registration_errors(final_transform, fixed_points, moving_points)\n",
    "\n",
    "print(optimizer_termination)\n",
    "print('After registration, errors in millimeters, mean(std): {:.2f}({:.2f}), max: {:.2f}'.format(final_errors_mean, final_errors_std, final_errors_max))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare the error distributions\n",
    "\n",
    "To fairly compare the two registration above we look at their running times (see results above) and their \n",
    "error distributions (plotted below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.hist(original_resolution_errors, bins=20, alpha=0.5, label='original resolution', color='blue')\n",
    "plt.hist(resampled_resolution_errors, bins=20, alpha=0.5, label='higher resolution', color='green')\n",
    "plt.legend()\n",
    "plt.title('TRE histogram');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "\n",
    "It appears that the memory-time trade-off works in our favor, but is this always the case? Well, you will have to answer that for yourself.\n",
    "\n",
    "Some immediate things you can try:\n",
    "* Change the interpolation method for the \"original resolution\" registration to nearest neighbor.\n",
    "* Change the resolution of the resampled image - will a higher resolution always result in faster running times?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
